function t = blockconv3d(u, Hx, Hy, Hz)
    t = real(ifftn(Hx.*fftn(u(:,:,:,1), size(Hx)) + ...
                   Hy.*fftn(u(:,:,:,2), size(Hy)) + ...
                   Hz.*fftn(u(:,:,:,3), size(Hz))));
end